PARTE 1 — Introducción

Cuarta esquina o fourth-corner

El problema de la cuarta esquina (fourth-corner) trata de relacionar los rasgos de los organismos con el entorno en el que estos se encuentran. Es una aproximación más holística al modelado de distribución de especies, que clásicamente se realiza teniendo en cuenta la relación directa de las especies con el ambiente. ¿Por qué recibe el nombre de la cuarta esquina? Se debe a que involucra tres matrices para calcular una cuarta, como vemos en el siguiente esquema:

Esquema modificado de Brown et al. (2014).
Esquema modificado de Brown et al. (2014).
  • La matriz L contiene información de especies por localidad
  • La matriz R contiene información del ambiente de dichas localidades
  • La matriz Q contiene información de los rasgos de las especies
  • La matriz D, entonces, contendría información de la relación rasgo-ambiente

Y, ¿por qué decimos que es un problema? Esto se debe a que no es nada sencillo calcular dicha matriz. En los últimos 30 años ha habido dos propuestas principales:

  • Análisis RLQ (Dolédec et al. 1996): consiste en un enfoque exploratorio que construye ordenaciones conjuntas de variables ambientales y rasgos mediante descomposición de valores singulares, proporcionando una visión cualitativa de asociaciones generales entre rasgos y ambiente (panel a de la figura a continuación)
  • Pruebas de hipótesis (Legendre et al. 1997): consiste en un enfoque de álgebra matricial en el que se realizan pruebas de permutación sobre las matrices R, L y Q para determinar una matriz D que resume qué asociaciones entre variables ambientales y rasgos son significativas, sin proporcionar detalles sobre la intensidad de estas interacciones (panel b de la figura a continuación)
Aproximaciones al problema de la cuarta esquina. Modificado de Brown et al. (2014).
Aproximaciones al problema de la cuarta esquina. Modificado de Brown et al. (2014).

Aunque ambas soluciones cumplen diferentes propósitos, no proporcionan una solución cuantitativa para la cuarta esquina y, por tanto, no nos permiten entender la magnitud de las interacciones entre rasgos y ambiente. En los pasos siguientes, veremos cómo podemos hacer esto utilizando modelos lineales generalizados con las tres matrices, siguiendo la solución propuesta por Brown et al. (2014).

Proceso de análisis (workflow)

Librerías de R

En primer lugar, cargaremos las librerías de R que vamos a utilizar a lo largo de nuestros análisis. Las librerías o paquetes son recopilaciones de funciones, generalmente sobre un temática específica, que nos sirven para realizar procesos en R, como manejo o análisis de datos, de una forma estandarizada y rápida. Aunque estos paquetes suelen estar escritos por personas que se dedican a la programación o a la ciencia, R es una comunidad de software libre, así que cualquier persona que tenga una metodología y quiera compartirla con la comunidad de usuarios de R puede crear su propio paquete. Los paquetes de R pueden instalarse directamente a través del CRAN utilizando la función install_packages() con el nombre del paquete como argumento de la función, entre comillas (por ejemplo: install.packages("tidyverse")). En ocasiones, también pueden instalarse directamente desde las cuentas de GitHub de los usuarios. Los paquetes solo necesitan instalarse una vez en un ordenador (aunque conviene estar informado sobre posibles actualizaciones). Sin embargo, cada vez que iniciamos una nueva sesión de R hemos de cargar esos paquetes. Esto se hace usando la función library(), que es la forma de decirle a nuestro programa que queremos utilizar las funciones que el paquete contiene.

A continuación, pues, instalamos las librerías a utilizar:

install.packages("tidyverse") # manipulación, transformación y dataviz
install.packages("mvabund") # manupilación y análisis de datos multivariantes

Y las cargamos:

# fíjate en que ya no hacen falta las comillas para cargar los paquetes
library(tidyverse)
library(mvabund)

Estos paquetes son:

  • tidyverse (Wickham et al. 2019), que es en realidad una colección de paquetes de programación funcional que sirven para la manipulación de datos de forma más eficiente, utilizando tuberías (pipes) para enlazar elementos. La pipe es %>%, que veremos mucho en el código en este documento. Aunque parezca poco tostón escribir %, > y %, podéis insertarla automáticamente presionando las teclas shift + control + m en Windows (o shift + command + m en Mac).

  • mvabund (Wang et al. 2012), que es un paquete para modelizar datos de abundancia en ecología de comunidades y que contiene las funciones que utilizaremos en el análisis de fourth-corner.

Datos dummy

Vamos a crear matrices con datos ficticios (lo que se llama en inglés datos dummy) para aprender cómo se utiliza esta técnica. Lo primero que haremos será seleccionar cuántas especies, cuántos rasgos funcionales, cuántas variables ambientales y cuántos sitios o localidades queremos:

species <- 20 # número de especies para nuestros datos dummy
sites <- 10 # número de sitios o localidades para nuestros datos dummy
traits <- 15 # número de rasgos funcionales para nuestros datos dummy
env_var <- 15 # número de variables ambientales para nuestros datos dummy

Y, a continuación, vamos a sembrar una semilla (seed). ¿Qué quiere decir esto? La semilla es un número natural que transforma lo estocástico en determinista. Esto es, al trabajar con funciones que generan datos aleatorios (por ejemplo, runif(n), siendo n el número de elementos a generar), si no hemos plantado una semilla, cada vez que corramos la función los resultados generados cambiarán (¡probad a correr runif(10) sin semilla varias veces!). Sin embargo, cuando plantamos una semilla usando la función set.seed(), el resultado generado será siempre el mismo, y sólo cambiará si cambiamos la semilla. Aquí utilizaremos el número 5558 como nuestra semilla por defecto.

set.seed(5558)

Finalmente, vamos a crear las matrices necesarias para los análisis. Tenemos en mente que, aunque en este ejercicio estamos generando datos ficticios, en la práctica estas matrices comprenderían datos reales de nuestro sistema de estudio. La cuestión es que estos datos deben estar organizados en las siguientes matrices con la siguiente estructura:

Matriz L

Matriz de abundancia de especies por localidad (L), con 20 especies en columnas y 10 localidades en filas.

# crear una base de datos de combinaciones de especie y sitio
L <- expand.grid(species = paste0("sp", 1:species),
                 sites = paste0("site", 1:sites))

# añadir una columna con los datos de abundancia
L$value <- rpois(nrow(L), 30)

# incorporar los datos a una matriz
L <- L %>%
  spread(key = species, value = value) %>%
  column_to_rownames(var = "sites") %>%
  as.data.frame()

# observar las primeras filas de dicha matriz
head(L)

Matriz Q

Matriz con datos de rasgos funcionales por especie (Q), con 15 rasgos en filas y 20 especies en columnas.

# crear una base de datos de combinaciones de rasgos y especies
Q <- expand.grid(species = paste0("sp", 1:species),
                 traits = paste0("t", 1:traits))

# añadir una columna con los datos cuantitativos
Q$value <- round(runif(n = nrow(Q), 0, 1), 3)

# incorporar los datos a una matriz
Q <- Q %>%
  spread(key = species, value = value) %>%
  column_to_rownames(var = "traits") %>%
  as.data.frame()

# observar las primeras filas de dicha matriz
head(Q)

2.2. Matriz R

Matriz con datos de variables ambientales por localidad (R), con 15 variables ambientales y 10 localidades.

# crear una base de datos de comb. de variables ambientales y sitios
R <- expand.grid(environ = paste0("env", 1:env_var),
                 sites = paste0("site", 1:sites))

# Add a column with random quantitative data from 0 to 1
R$value <- round(runif(n = nrow(R), 0, 1), 3)

# incorporar los datos a una matriz
R <- R %>%
  spread(key = environ, value = value) %>%
  column_to_rownames(var = "sites") %>%
  as.data.frame()

# observar las primeras filas de dicha matriz
head(R)

Exploración de los datos

El objetivo de esta sesión no es aprender R, ni mucho menos dominar la visualización de datos (que requiere mucha práctica), pero sí me gustaría que os lleváseis una idea de que estas herramientas están disponibles y no son tan innaccesibles como pueden parecer.

Lo primero que vamos a hacer es observar cómo es la distribución de la abundancia de organismos en las diferentes poblaciones. Para ello, vamos a utilizar el paquete ggplot2, perteneciente al tidyverse, del que ya hablamos anteriormente. La función ggplot() es una de las más poderosas del lenguaje R para crear visualizaciones de datos. Funciona por capas: primero, ggplot crea un espacio o lienzo definido por ciertas variables de una base de datos. Dicho lienzo puede ser completado con la adición de distintas funciones de representación de los datos, como geom_point() o geom_line() que añaden una nube de puntos o una línea, respectivamente. Asimismo, funciones cómo theme() se encargan de cuidar los detalles estéticos del gráfico, como el tamaño de la letra de los rótulos, el color de fondo o la presencia o ausencia de ejes. Para más comodidad, en lugar de añadir la estética para cada gráfico, podemos utilizar la función theme_set() para definir una estética global que se aplica a todos los gráficos.

# estética general de los plots
theme_set(theme_minimal() +
            theme(axis.title = element_text(size = 15),
                  axis.line = element_line(color = "black", linewidth = 0.5),
                  panel.grid.major = element_blank(), 
                  panel.grid.minor = element_blank(),
                  axis.text = element_text(color = "black", size = 12),
                  strip.text.x = element_text(size = 12),
                  axis.ticks = element_line(color = "black")))
# plot de abundancia por sitio
L %>%
  
  # pasamos los nombres de fila a una columna
  rownames_to_column(var = "site") %>% 
  
  # transformamos la matrix para poder hacer el plot
  pivot_longer(cols = -site,
               names_to = "species",
               values_to = "abundance") %>% 
  
  # añadimos el plot
  ggplot(aes(x = abundance,
             y = reorder(species, abundance, FUN = max),
             color = site)) +
  geom_point() + # gráfico de puntos
  labs(x = "Abundancia",
       y = "Organismo",
       color = "Sitio") +
  theme(text = element_text(size = 18))

Aproximaciones a la distribución de especies

Queremos entender cómo se distribuyen los organismos en términos de abundancia. A continuación, vamos a ir construyendo nuestro entendimiento de cómo esto sucede, desde las preguntas y análisis más sencillos.

¿Efecto general del sitio?

Lo primero de todo: hemos muestreado diez localidades en cinco espacios naturales diferentes. ¿Está la abundancia de los diferentes organismos determinada por el lugar de muestreo?

El paquete mvabund nos permite ajustar modelos para explicar la abundancia de los organismos. Para esto, es necesaria una pequeña transformación de la estructure de la matriz de datos L, que se realiza con la función homónima mvabund(), integrada en el paquete. Tras esto, somos capaces de utilizar otra función de mvabund, manyglm(), para explicar la abundancia de la matriz en función del lugar de muestreo, explícito como rownames(L) o, lo que es lo mismo, los nombres de las filas de la matriz L. Podemos especificar la distribución de los datos para el modelo utilizando el argumento family; en este caso, lo que mejor se ajusta es una distribución poisson, pero podéis probar a cambiar esta distribución (por ejemplo, a family = "negative_binomial") y ver qué ocurre con la distribución de los residuos.

# creamos un objeto 'mvabund' a partir de la abundancia
spp <- mvabund(L)

# hacemos un modelo para explicar la abundancia en función del sitio
mod <- manyglm(spp ~ rownames(L), family = "poisson") # distr. poisson

# representamos los residuos del modelo
plot(mod)

Luego, simplemente lo podemos resolver con un análisis de la varianza:

# análisis de la varianza para ver si el sitio explica la abundancia
(anov <- anova(mod, nBoot = 999))
## Time elapsed: 0 hr 0 min 1 sec
## Analysis of Deviance Table
## 
## Model: spp ~ rownames(L)
## 
## Multivariate test:
##             Res.Df Df.diff   Dev Pr(>Dev)
## (Intercept)      9                       
## rownames(L)      0       9 156.5    0.201
## Arguments:
##  Test statistics calculated assuming uncorrelated response (for faster computation) 
##  P-value calculated using 999 iterations via PIT-trap resampling.

Como podemos ver, el estadístico nos indica que p = 0.201 > 0.05 y, por tanto, podemos decir que el sitio no tiene un efecto significativo a la hora de explicar las abundancias de los organismos.

Explicación a través del ambiente

En segundo lugar, podemos realizar la aproximación más clásica en modelos de distribución de especies: intentar entender qué variables climáticas explican la abundancia de los diferentes organismos. Ciertos modelos más o menos complejos incluso permiten hacer predicciones sobre la distribución potencial de los organismos a través de modelización ambiental.

En nuestro caso, podemos utilizar una función integrada en el paquete mvabund llamada traitglm(). Esta función aplica modelos lineales generalizados (GLMs) para explicar una matriz de datos en función de otra (u otras) matrices de datos de estructura compatible. Por ejemplo, no podremos explicar una matriz de 10x20 en función de una matriz de 5x15, pero sí una de 10x20 en función de una de 10x15.

# explicamos la matriz de abundancia en función del ambiente
sdm_env <- traitglm(L, R) # modelo lineal generalizado (glm)
## No traits matrix entered, so will fit SDMs with different env response for each spp
# transformar y dataviz
sdm_env$fourth.corner %>%
  as_tibble(rownames = "species") %>%
  pivot_longer(cols = -species, names_to = "env_var", values_to = "value") %>% 
  mutate(species = as.factor(str_extract(species, "sp\\d+")),
         value = round(value, 5),
         env_var = forcats::fct_inorder(env_var),
         species = forcats::fct_reorder(species, as.numeric(str_extract(species, "\\d+")))) %>% 

  # mapa de calor
  ggplot(aes(x = env_var, y = fct_rev(species), fill = value)) +
  geom_raster() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red",
    midpoint = 0) +
  labs(x = "Variables ambientales",
       y = "Especies",
       fill = "Coeficiente") +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))

A pesar de que esta aproximación nos aporta gran información sobre la relación de las especies con su ambiente, no deja de ser una forma de analizar qué ocurre con las especies de forma individual. Esto no nos permite alcanzar el nivel de abstracción que deseamos para resolver a nuestra pregunta; para ello, necesitamos incluir los datos de los rasgos funcionales de las especies.

Rasgos funcionales y la matriz fourth-corner

Esto lo podemos hacer con la misma función, traitglm(), pero añadiendo la matriz Q. En este caso, necesitamos introducir la transpuesta (t()) de Q porque sino el modelo no puede computar matrices de diferente tamaño y nos da el siguiente error: number of rows of matrices must match (see arg 2).

# añadimos la matriz de rasgos funcionales
f.corner <- traitglm(L, R, t(Q)) # glm de nuevo

# observemos la matriz fourth-corner
f.corner$fourth.corner
##             env1         env2         env3         env4         env5
## t1   0.005896477 -0.030156079  0.032728336 -0.022556718 -0.026917732
## t2  -0.015514887  0.013106324 -0.012049221 -0.039692052 -0.007347624
## t3   0.020185284  0.024309594  0.049111168  0.001369728 -0.046107386
## t4  -0.004948459  0.012092828 -0.030603234 -0.014503910  0.024181498
## t5  -0.020400040  0.005302025  0.020634272 -0.017760382 -0.009821665
## t6  -0.042718373  0.017132002 -0.056464978 -0.041663403  0.040201806
## t7   0.016447268 -0.032752837  0.066954612  0.021185391 -0.057901879
## t8   0.006364846 -0.025288544  0.051488167  0.018310605 -0.029164103
## t9   0.014388225 -0.047493117  0.024863594  0.010236594  0.014282519
## t10 -0.013240487  0.036835328 -0.015586492 -0.008098365 -0.008877121
## t11 -0.014959676 -0.013503966  0.008225036  0.007433749 -0.012406317
## t12  0.003445801 -0.036153152  0.051730749  0.009830930 -0.019343220
## t13 -0.019283181 -0.024839259  0.048699685 -0.016277801 -0.036393665
## t14 -0.024046745  0.050593395 -0.028186839 -0.009237108  0.013941703
## t15 -0.004714214 -0.001967574 -0.007937043 -0.026548242  0.027757722
##             env6         env7         env8         env9        env10
## t1  -0.019726434 -0.029070770  0.015378237 -0.033377981 -0.007832847
## t2   0.001187777 -0.021322501  0.035382373 -0.026892351  0.015762487
## t3   0.028208979 -0.008266815 -0.008361714  0.011775754  0.032692182
## t4   0.004730038  0.018637453  0.007642035 -0.002224654 -0.015310601
## t5  -0.015122491  0.024281807  0.001866518 -0.009335415 -0.019099460
## t6  -0.020126082  0.072155957  0.013309618  0.008538976 -0.026845870
## t7  -0.025845623 -0.002285457 -0.034460053  0.008901942 -0.005117371
## t8  -0.010154626 -0.037706166 -0.010054318 -0.012561046  0.007880418
## t9   0.010102239 -0.080896695  0.030564598 -0.030114349  0.037166281
## t10  0.002662576  0.038930905 -0.002267904  0.016424265 -0.003331159
## t11 -0.033098305  0.024245015 -0.013565228  0.009772219 -0.013308426
## t12 -0.012807564 -0.041765937  0.031867364 -0.015526612  0.037339151
## t13 -0.040957657  0.010185582 -0.017816583 -0.018358099 -0.031138063
## t14  0.001621939  0.028463857  0.031573752  0.001838934  0.005727427
## t15  0.010479725  0.005533823  0.012901708 -0.014480981 -0.010271220
##            env11        env12        env13        env14        env15
## t1   0.028755140  0.027739712 -0.001766501 -0.008262744 -0.011906722
## t2   0.018141164  0.050584237  0.034553432  0.011548849 -0.013527046
## t3   0.001882132  0.004668550  0.025864721 -0.003212201  0.030810105
## t4  -0.006084416  0.009707665 -0.011986105 -0.005002110 -0.042549437
## t5   0.004529213  0.019836164  0.011786953 -0.006707977 -0.024753065
## t6   0.006919360  0.036652412  0.013557160  0.022501647 -0.060949812
## t7   0.015657700 -0.018618540 -0.001088288  0.003257411  0.038172591
## t8   0.011025947 -0.015608570  0.010151306 -0.001234242  0.053281915
## t9   0.014327452 -0.009588088 -0.006568942 -0.020472294  0.065077764
## t10 -0.015113117  0.015529544  0.012788965  0.003663808 -0.032685636
## t11  0.013128807 -0.003697762  0.009224238  0.023593957  0.012053497
## t12  0.026922428  0.007788756  0.008441355 -0.009807739  0.055847632
## t13  0.026987546  0.014832604  0.021859639  0.013164277  0.003823548
## t14  0.001589585  0.027440968  0.024061261  0.018156810 -0.030603694
## t15  0.006159755  0.014595450 -0.005079241 -0.012195386 -0.026284344

Y, de esta forma, somos capaces de calcular las interacciones entre los rasgos funcionales de los organismos y el ambiente en el que ocurren. Vamos a representar estos coeficientes gráficamente, no sin antes hacer una rápida comprobación de los residuos del modelo:

# comprobamos las asunciones del modelo
plot(f.corner)

# finalmente, podemos representar el mapa de calor
f.corner$fourth.corner %>%
  as_tibble(rownames = "traits") %>%
  pivot_longer(cols = -traits, names_to = "env_var", values_to = "value") %>% 
  mutate(value = round(value, 5),
         env_var = forcats::fct_inorder(env_var),
         traits = forcats::fct_reorder(traits, as.numeric(str_extract(traits, "\\d+")))) %>% 
  
  # mapa de calor
  ggplot(aes(x = env_var, y = fct_rev(traits), fill = value)) +
  geom_raster() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red",
    midpoint = 0) +
  labs(x = "Variables ambientales",
       y = "Rasgos funcionales",
       fill = "Coeficiente") +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))

¿Qué observamos en esta matriz? Se ve que todas las posibles interacciones han sido otorgadas un determinado valor. ¿Es esto realista? ¿Son todas las interacciones entre rasgos funcionales y variables ambientales importantes?

Regularización lasso

Lasso (least square shrinkage and selection operator) (Tibshirani 1996) es una forma de regresión que añade una penalización a los coeficientes de mínimos cuadrados, que es el cálculo que se utiliza para encontrar modelos de regresión por comparación de los valores predichos con los observados. Esta penalización puede modificar los coeficientes de las variables del modelo y está controlada por un parámetro, \(\lambda\), que puede tomar valores desde 1 hasta infinito positivo. A medida que \(\lambda\) aumenta, la pendiente de la variable en el modelo mengua. La grandísima ventaja de esto es que lasso es capaz de llevar a cero la pendiente de aquellas variables del modelo que no sean relevantes para explicar los datos, dejando solo aquellas que sí lo sean. Esto nos permite simplificar los modelos y, en general, nuestro trabajo y toma de decisiones. Como recurso, os dejo este vídeo sobre regularización lasso que, aunque en un tono no demasiado serio, consigue explicar bien qué utilidad tiene y cómo funciona. Es ese vídeo oiréis hablar también de la ridge regression, que es muy similar a la regresión lasso pero no es capaz de llevar los coeficientes de las variables poco explicativas a cero.

Lasso está disponible para la aproximación fourth-corner a través del algoritmo glm1path de la función traitglm(), como se puede ver a continuación. Sin embargo, lasso no es un método exclusivo de fourth-corner, sino que se puede aplicar a cualquier regresión que necesite selección de variables, como por ejemplo en el análisis de big data.

# añadimos lasso a través del argumento method = "glm1path"
lasso <- traitglm(L, R, t(Q), method = "glm1path")

# observemos la matriz fourth-corner con lasso
lasso$fourth.corner
##     env1        env2          env3 env4 env5 env6        env7 env8 env9 env10
## t1     0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
## t2     0 0.000000000 -0.0126001270    0    0    0 0.000000000    0    0     0
## t3     0 0.000000000  0.0007411432    0    0    0 0.000000000    0    0     0
## t4     0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
## t5     0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
## t6     0 0.000000000  0.0000000000    0    0    0 0.009344633    0    0     0
## t7     0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
## t8     0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
## t9     0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
## t10    0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
## t11    0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
## t12    0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
## t13    0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
## t14    0 0.006611215  0.0000000000    0    0    0 0.000000000    0    0     0
## t15    0 0.000000000  0.0000000000    0    0    0 0.000000000    0    0     0
##     env11 env12 env13 env14       env15
## t1      0     0     0     0 0.000000000
## t2      0     0     0     0 0.000000000
## t3      0     0     0     0 0.000000000
## t4      0     0     0     0 0.000000000
## t5      0     0     0     0 0.000000000
## t6      0     0     0     0 0.005031213
## t7      0     0     0     0 0.000000000
## t8      0     0     0     0 0.000000000
## t9      0     0     0     0 0.000000000
## t10     0     0     0     0 0.000000000
## t11     0     0     0     0 0.000000000
## t12     0     0     0     0 0.000000000
## t13     0     0     0     0 0.000000000
## t14     0     0     0     0 0.000000000
## t15     0     0     0     0 0.000000000
# y volvemos a representar el mapa de calor
lasso$fourth.corner %>%
  as_tibble(rownames = "traits") %>%
  pivot_longer(cols = -traits, names_to = "env_var", values_to = "value") %>% 
  mutate(value = round(value, 5),
         env_var = forcats::fct_inorder(env_var),
         traits = forcats::fct_reorder(traits, as.numeric(str_extract(traits, "\\d+")))) %>% 
  
  # mapa de calor
  ggplot(aes(x = env_var, y = fct_rev(traits), fill = value)) +
  geom_raster() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red",
    midpoint = 0) +
  labs(x = "Variables ambientales",
       y = "Rasgos funcionales",
       fill = "Coeficiente") +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))

Ejercicio — Parte 1

La semilla que sembramos al principio, 5558, está seleccionada para generar unos datos que producen estos resultados concretos. ¿Qué ocurre si cambiamos la semilla? Probad a rehacer los análisis con semillas que tomen los siguientes valores: 7120, 9937, 1840, 5647 y 216.

PARTE 2 — Aplicación a sistemas fluviales

Base de datos

En esta segunda parte vamos a dejar de lado los datos dummy y, en su lugar, vamos a analizar una base de datos real de macroinvertebrados de agua dulce. Se trata de 76 familias de organismos, desde efímeras hasta gammáridos. La lista completa se puede encontrar usando el siguiente comando: colnames(L), siempre después de haber cargado la matriz L correspondiente (siguiente paso). Asimismo, los datos proceden de cinco diferentes espacios naturales de los Montes Vascos, a saber: AIZK = Aizkorri (Gipuzkoa), ARAL = Aralar (Gipuzkoa), ARTI = Artikutza (Nafarroa), GORB = Gorbeia (Bizkaia/Araba), e IZKI = Izki (Araba).

Datos de abundancia por localidad (matriz L)

Cargamos la matriz L con el siguiente comando:

L <- read.table("data/L_family.txt")
head(L)

Datos de rasgos por familia (matriz Q)

Cargamos ahora la matriz Q de rasgos funcionales por familia:

Q <- read.table("data/Q_family.txt")
head(Q)

En la tabla siguiente se muestra la descripción de los rasgos funcionales de los organismos. A tener en cuenta: los valores de los rasgos se corresponden con los valores mostrados como columnas. Esto quiere decir que, por ejemplo, si la variable 1 (var1) tiene un valor de 1.33, ese organismo puede alcanzar un tamaño máximo de entre 0.25cm y 5cm, pero es más probable que su tamaño sea de 0.25cm o menor (porque 1.33 está más cerca de 1 que de 2).

Rasgos funcionales de los organismos:
Valor del rasgo
1 2 3 4 5 6 7
Tamaño potencial máximo del organismo
var1 X0.25cm X0.250.5cm X0.51cm X12cm X24cm X48cm X8cm
Distribución longitudinal
var13 crenon epirithron metarithron hyporithron epipotamon metapotamon estuary
Altitud en la que se encuentra
var14 lowlands piedmont.level alpine.level
Velocidad de corriente (preferencia)
var17 null slow.less.25.cms medium.25.50.cmlesss fast.more.50.cms
Status trófico (preferencia)
var18 oligotrophic mesotrophic eutrophic
Saprobicidad (preferencia)
var21 xenosaprobic oligosaprobic b.mesosaprobic a.mesosaprobic polysaprobic
pH del medio (preferencia)
var22 less.4 more.4.4.5 more.4.5.5 more.5.5.5 more.5.5.6 more.6

 

Datos de variables ambientales por localidad (matriz R)

Y, finalmente, cargamos la matriz R con las variables ambientales:

R <- read.table("data/R_family.txt")
head(R)

En la siguiente tabla se proporciona información sobre el significado de las variables ambientales medidas en los espacios naturales muestreados.

Explicación de las variables ambientales utilizadas:
Variable Descripción
pH Valores de pH en el medio
Conductivity_.µS.cm. Conductividad (microS/cm)
O2_sat Saturación de oxígeno (%)
O2_conc Concentración de oxígeno (mg/L)
Cl_.mg.L._avg Concentración de cloruro (mg/L)
SO4_.mg.L._avg Concentración de sulfato (mg/L)
NO3_.mg.L._avg Concentración de nitrato (mg/L)
N.NH4_.mg.L._avg Concentración de amonio (mg/L)
DIN.N.mg.L. Concentración de nitrógeno (mg/L)
P.PO4_.µg.L. Concentración de fosfato (mg/L)
tree_cover_. Cobertura arbórea (%)
BARE_SOIL_._19_23 Cobertura de suelo descubierto (%)
GRASSLAND_._19_23 Cobertura de praderas (%)
SHRUBLAND_._19_23 Cobertura arbustiva (%)
NAT_FOREST_._19_23 Cobertura de bosque natural (%)
PLANT_FOREST_._19_23 Cobertura de plantaciones (%)
Catchment_area_km2 Área de muestreo (km2)

 

Análisis de datos

El objetivo es calcular las interacciones significativas entre los rasgos funcionales y las variables ambientales que nos ayuden a explicar la distribución de la abundancia de los diferentes organismos. Para ello, comenzaremos visualizando los datos y, después, utilizaremos la función traitglm() del paquete mvabund tal y como vimos anteriormente: primero, sin penalización lasso y, segundo, con ella para seleccionar interacciones.

Exploración de los datos

Vamos a hacer un par de gráficos que nos ayuden a entender mejor qué datos tenemos.

# plot de abundancia por sitio
L %>%
  
  # rownames (sites) to column to plot
  rownames_to_column(var = "site") %>% 
  
  # transformamos la matrix para poder hacer el plot
  pivot_longer(cols = -site,
               names_to = "species",
               values_to = "abundance") %>% 
  
  # añadimos el plot
  ggplot(aes(x = abundance,
             y = reorder(species, abundance, FUN = max),
             color = site)) +
  geom_point() + # gráfico de puntos
  labs(x = "Abundancia",
       y = "Familias",
       color = "Sitio") +
  theme(text = element_text(size = 18))

¿Se parece esta distribución de las abundancias a la que teníamos anteriormente para los datos dummy?

Interacción rasgos-ambiente, con y sin penalización lasso

Vamos ahora directamente a resolver la pregunta que nos interesa, para lo cual utilizaremos las tres matrices como argumentos de traitglm(), y representamos gráficamente la matriz de interacciones seleccionada. Lo hacemos para ambos modelos, sin lasso y con lasso (recordamos: method = glm1path) y los representamos.

Primero sin lasso:

f.corner <- traitglm(L, R, Q)

f.corner$fourth.corner %>%
  as_tibble(rownames = "traits") %>%
  pivot_longer(cols = -traits, names_to = "env_var", values_to = "value") %>% 
  mutate(value = round(value, 5),
         env_var = factor(env_var, levels = sort(unique(env_var))),
         traits = factor(traits, levels = sort(unique(traits), decreasing = TRUE))) %>% 
  
  # mapa de calor
  ggplot(aes(x = env_var, y = traits, fill = value)) +
  geom_raster() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red",
                       midpoint = 0) +
  labs(title = "Modelo independiente del espacio natural",
       x = "Var. ambientales",
       y = "Rasgos func.",
       fill = "Coef.") +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))

Y después con lasso:

lasso <- traitglm(L, R, Q, method = "glm1path")

lasso$fourth.corner %>%
  as_tibble(rownames = "traits") %>%
  pivot_longer(cols = -traits, names_to = "env_var", values_to = "value") %>% 
  mutate(value = round(value, 5),
         env_var = factor(env_var, levels = sort(unique(env_var))),
         traits = factor(traits, levels = sort(unique(traits), decreasing = TRUE))) %>% 
  
  # mapa de calor
  ggplot(aes(x = env_var, y = traits, fill = value)) +
  geom_raster() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red",
                       midpoint = 0) +
  labs(title = "Modelo independiente del espacio natural",
       x = "Var. ambientales",
       y = "Rasgos func.",
       fill = "Coef.") +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))

¿Qué está pasando aquí? ¿Por qué parece que hay menos interacciones seleccionadas para el gráfico sin lasso (izquierda) que para el gráfico con lasso (derecha)? Vamos a ver las matrices de interacciones:

print(f.corner$fourth.corner)
##                pH Conductivity_.µS.cm.      O2_sat     O2_conc Cl_.mg.L._avg
## var1   0.01315904           0.08485679 -0.04385779 -0.09429455   0.159752592
## var13  0.03724744           0.07064226 -0.12026136  0.10177320  -0.130719892
## var14  0.01182520          -0.20460362 -0.06799693  0.05762894  -0.068914605
## var17 -0.22817063           0.03082527  0.36298936 -0.05095678  -0.308052877
## var18 -0.04010953           0.09273085  0.12981674 -0.03694424  -0.066904964
## var21 -0.15739570          -0.01029542  0.05778800 -0.14558510   0.037993393
## var22 -0.14464822           0.08830754  0.09045460 -0.02330884  -0.003115431
##       SO4_.mg.L._avg NO3_.mg.L._avg N.NH4_.mg.L._avg DIN.N.mg.L. P.PO4_.µg.L.
## var1     -0.13424086    -0.07621326      0.038806233 -0.07479167  -0.12263507
## var13    -0.01803699     0.13607971      0.210911512  0.14646150   0.01439262
## var14     0.10573982    -0.03919476      0.026474914 -0.03816324   0.02472374
## var17    -0.05487787     0.08396887      0.020640364  0.08532485  -0.14851299
## var18    -0.12750873    -0.08065156      0.001747011 -0.08096013  -0.11096515
## var21     0.05871555    -0.03579921     -0.045194207 -0.03805644   0.04864642
## var22    -0.08757983     0.07776019     -0.260557410  0.06612103   0.10116822
##       tree_cover_. BARE_SOIL_._19_23 GRASSLAND_._19_23 SHRUBLAND_._19_23
## var1   0.079114951       -0.11543830        -0.1971582        -0.2500568
## var13  0.012928469        0.29542957         0.9264349         1.7104471
## var14 -0.005695622       -0.06749943        -0.3468541        -0.8438421
## var17  0.041052274       -0.01778968         0.2630383         0.9119774
## var18  0.064415597        0.02771697         0.3715773         0.7648960
## var21 -0.060071252       -0.18598408        -0.5006729        -1.3722988
## var22 -0.122324576       -0.06634006        -0.1087394        -0.4405074
##       NAT_FOREST_._19_23 PLANT_FOREST_._19_23 Catchment_area_km2
## var1          -0.4186024           -0.1426231        -0.01484885
## var13          2.4691483            1.7927658         0.06448584
## var14         -1.1708913           -0.8676686        -0.15495297
## var17          1.4136504            0.9392763         0.22122586
## var18          1.3579632            1.0860101         0.05005877
## var21         -1.8517342           -1.3649009         0.06347529
## var22         -0.7797123           -0.4636749         0.01265808
print(lasso$fourth.corner)
##                pH Conductivity_.µS.cm.      O2_sat      O2_conc Cl_.mg.L._avg
## var1   0.00000000           0.04532127 -0.04740023 -0.040929013    0.13017623
## var13  0.00000000           0.00000000 -0.03418713  0.000000000   -0.02853290
## var14  0.00000000          -0.07377945 -0.02252036  0.000000000   -0.07853568
## var17 -0.06788648          -0.13447846  0.15216154  0.015706817   -0.04396393
## var18  0.00000000           0.00000000  0.00000000  0.000000000    0.00000000
## var21 -0.01216116           0.02524483  0.00000000 -0.054184873    0.02458062
## var22 -0.02059854           0.00000000  0.06917813  0.003542237    0.00000000
##       SO4_.mg.L._avg NO3_.mg.L._avg N.NH4_.mg.L._avg DIN.N.mg.L. P.PO4_.µg.L.
## var1               0    -0.09755478      0.016320134  0.00000000 -0.069472661
## var13              0     0.00000000      0.076288129  0.07944128  0.020166768
## var14              0     0.00000000      0.001604626  0.00000000  0.000000000
## var17              0     0.01974271      0.000000000  0.00000000 -0.004395314
## var18              0     0.00000000      0.000000000 -0.10166980  0.000000000
## var21              0     0.00000000      0.000000000  0.00000000  0.000000000
## var22              0     0.06100831     -0.153173000  0.00000000  0.022676177
##       tree_cover_. BARE_SOIL_._19_23 GRASSLAND_._19_23 SHRUBLAND_._19_23
## var1    0.03139716       -0.03484952       -0.08209893       0.000000000
## var13   0.00000000        0.00000000        0.16596512       0.000000000
## var14   0.01500141        0.00000000        0.00000000       0.000000000
## var17   0.00000000       -0.05502395       -0.03508282       0.002630367
## var18   0.00000000       -0.01780711        0.00000000      -0.056330545
## var21  -0.02603755        0.00000000        0.00000000      -0.069400967
## var22  -0.03645761        0.00000000        0.03377143       0.000000000
##       NAT_FOREST_._19_23 PLANT_FOREST_._19_23 Catchment_area_km2
## var1        -0.021488948           0.02589767       -0.009847158
## var13        0.000000000           0.00000000        0.093095220
## var14        0.000000000           0.00000000       -0.125846106
## var17        0.000000000          -0.01222346        0.136464856
## var18        0.001169351           0.00000000        0.000000000
## var21        0.000000000           0.00000000        0.000000000
## var22       -0.037800510           0.00000000        0.034095625

Como vemos, hemos de ir con ojo con las representaciones visuales, ya que una escala de colores puede darnos una idea errónea de los resultados.

Ejercicio — Parte 2

El objetivo de este ejercicio final es resolver la siguiente pregunta: ¿qué ocurre con las interacciones rasgo-ambiente si analizamos cada espacio natural por separado?

Pista: primero tenemos que seleccionar solamente las filas que correspondan al espacio natural que nos interese. Esto lo podemos hacer con la función grepl() del siguiente modo: X[grepl("CODE", rownames(X)), ], siendo X la matriz que corresponda y siendo CODE el código del espacio natural deseado. El resto del proceso es idéntico a lo explicado anteriormente: traitglm() con las tres matrices y el argumento para usar lasso.

Aquí los resultados:

Referencias

Brown, A.M., Warton, D.I., Andrew, N.R., Binns, M., Cassis, G. & Gibb, H. (2014). The fourth-corner solution using predictive models to understand how species traits interact with the environment. Methods in Ecology and Evolution, 5, 344–352.
Dolédec, S., Chessel, D., Braak, C.J.F. ter & Champely, S. (1996). Matching species traits to environmental variables: a new three-table ordination method. Environmental and Ecological Statistics, 3, 143–166.
Legendre, P., Galzin, R. & Harmelin-Vivien, M.L. (1997). Relating Behavior to Habitat: Solutions to Thefourth-Corner Problem. Ecology, 78, 547–562.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methodological), 58, 267–288.
Wang, Y., Naumann, U., Wright, S.T. & Warton, D.I. (2012). mvabund an R package for model-based analysis of multivariate abundance data. Methods in Ecology and Evolution, 3, 471–474.
Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L.D., François, R., et al. (2019). Welcome to the Tidyverse. Journal of Open Source Software, 4, 1686.

  1. Basque Center for Climate Change (BC3); ↩︎